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Differential Complexes and 
Numerical Stability 



Douglas N. Arnold* 

Abstract 

Differential complexes such as the de Rham complex have recently come 
to play an important role in the design and analysis of numerical methods for 
partial differential equations. The design of stable discretizations of systems 
of partial differential equations often hinges on capturing subtle aspects of 
the structure of the system in the discretization. In many cases the differen- 
tial geometric structure captured by a differential complex has proven to be a 
key element, and a discrete differential complex which is appropriately related 
to the original complex is essential. This new geometric viewpoint has pro- 
vided a unifying understanding of a variety of innovative numerical methods 
developed over recent decades and pointed the way to stable discretizations 
of problems for which none were previously known, and it appears likely to 
play an important role in attacking some currently intractable problems in 
numerical PDE. 
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1. Introduction 

During the twentieth century chain complexes, their exactness properties, and 
commutative diagrams involving them pervaded many branches of mathematics, 
most notably algebraic topology and differential geometry. Recently such homolog- 
ical techniques have come to play an important role in a branch of mathematics 
often thought quite distant from these, numerical analysis. Their most significant 
applications have been to the design and analysis of numerical methods for the 
solution of partial differential equations. 

Let us consider a general problem, such as a boundary value problem in partial 
differential equations, as an operator equation: given data / in some space Y find 
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the solution u in some space X to the problem Lu = f. A numerical method 
discretizes this problem through the construction of an operator Lh : Xh Yh 
and data fh S Yh and defines an approximate solution Uh G X^ by the equation 
LhUh = fh- Of course the numerical method is not likely to be of value unless it is 
consistent which means that and fh should be close to L and / in an appropriate 
sense. 

Before we speak of solving the original problem, numerically or otherwise, we 
should first confront the question of whether it is wcU-posed. That is, given / G Y, 
does a unique u £ X exist, and, if so, do small changes of / induce small changes 
in y? The analogous questions for the numerical method, whether given fh S Yh 
a unique Uh € Xh is determined by the discrete equation LhUh = fh, and whether 
small changes in fh induce small changes in Uh-, is the question of stability of the 
numerical method. A common paradigm, which can be formalized in many contexts 
of numerical analysis, is that a method which is consistent and stable is convergent. 

Well-posedness is a central issue in the theory of partial differential equations. 
Of course, we do not expect just any PDE problem to be well-posed. Well-posedness 
hinges on structure of the problem which may be elusive or delicate. Superficially 
small changes, for example to the sign of a coefficient or the type of boundary con- 
ditions, can certainly destroy well-posedness. The same is true for the stability of 
numerical methods: it often depends on subtle or elusive properties of the numerical 
scheme. Usually stability reflects some portion of the structure of the original prob- 
lem that is captured by the numerical scheme. However in many contexts it is not 
enough that the numerical scheme be close to the original problem in a quantitative 
sense for it to inherit stability. That is, it may well happen that a consistent method 
for a well-posed problem is unstable. In this paper we shall see several examples 
where the exactness properties of discrete differential complexes and their relation 
to differential complexes associated with the PDE are crucial tools in establishing 
the stability of numerical methods. In some cases the homological arguments have 
served to elucidate or validate methods that had been developed over the preceding 
decades. In others they have pointed the way to stable discretizations of problems 
for which none were previously known. They will very likely play a similar role in 
the eventual solution of some formidable open problems in numerical PDE, espe- 
cially for problems with significant geometric content, such as in numerical general 
relativity. As in other branches of mathematics, in numerical analysis differential 
complexes serve both to encode ken- structure concisely and to unify considerations 
from seemingly very different contexts. 

In this paper we shall discuss only finite element methods since, of the major 
classes of numerical methods for PDE, they are the most amenable to rigorous 
analysis, and have seen the greatest use of differential complexes. But complexes 
have recently arisen in the study of finite differences, finite volumes, and spectral 
methods as well. 



2. Finite element spaces 

A finite element space on a domain O is a function space defined piecewise 
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by a certain assembly procedure which we now recall; cf. For simplicity, here 
we shall restrict to spaces of piecewise polynomials with respect to a triangulation 
of an n-dimensional domain by n-simplices with n = 2 or 3 (so implicitly we are 
assuming that f2 C is polygonal or C M'^ is polyhedral). On each simplex T 
we require that there be given a function space of shape function Wt and a set of 
degrees of freedom, i.e., a set of linear functionals on Wt which form a basis for 
the dual space. Moreover, each degree of freedom is supposed to be associated with 
some subsimplex of some dimension, i.e., in three dimensions with a vertex, an edge, 
a face, or the tetrahedron itself. For a subsimplex which is shared by two simplices 
in the triangulation, we assume that the corresponding functionals arc in one-to- 
one-correspondence. Then the finite element space Wh is defined as those functions 
on SI whose restriction to each simplex T of the triangulation belongs to Wt and for 
which the corresponding degrees of freedom agree whenever a subsimplex is shared 
by two simplices. 

The simplest example is obtained by choosing Wt to be the constant functions 
and taking as the only degree of freedom on T the 0th order moment (/) i— > (I'i^) 
(which we associate with T itself). The resulting finite element space is simply 
the space of piecewise constant functions with respect to the given triangulation. 
Similarly we could choose Wt = ^i{T) (by Fp(T) we denote the space of polynomial 
functions on T of degree at most p), and take as degrees of freedom the moments 
of degrees and also those of degree 1, </> i— > jj, (j){x)xi dx. Again all the degrees of 
freedom are associated to T itself. This time the finite element space consists of all 
piecewise linear functions. Of course, the construction extends to higher degrees. 

A more common piecewise linear finite element space occurs if we again choose 
Wt — Pi (T) , but take as degrees of freedom the maps (p i-^ (t){v) , one associated 
to each vertex v. In this case the assembled finite element space consists of all 
continuous piecewise linear functions. More generally we can choose Wt — fp{T) 
for p > 1, and associate to each vertex the evaluation degrees of freedom just 
mentioned, to each edge the moments on the edge of degree at most p ~ 2, to each 
face the moments on the face of degree at most p — 3, and to each tetrahedron 
the moments of degree at most p — A. The resulting finite element space, called the 
Lagrange finite element of degree p, consists of all continuous piecewise polynomials 
of degree at most p. Figure |l| shows a mesh of a two dimensional domain and a 
typical function in the space of Lagrange finite elements of degree 2 with respect to 
this mesh. 

Mnemonic diagrams as in Figure ^ are often associated to finite element spaces, 
depicting a single element T and a marker for each degree of freedom. 

Next we describe some finite element spaces that can be used to approximate 
vector-valued functions. For brevity we limit the descriptions to the 3-dimensional 
case, but supply diagrams in both 2 and 3 dimensions. Of course we may simply take 
the Cartesian product of three copies of one of the previous spaces. For example, 
the element diagrams shown on the left of Figure ^ refer to continuous piecewise 
linear vector fields in two and three dimensions. More interesting spaces are the 
face elements and edge elements essentially conceived by Raviart and Thomas 
in two dimensions and by Nedelec in three dimensions. In the lowest order 
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Figure 1: A mesh marked with the locations of the degrees of 
freedom for Lagrange finite elements of degree 2 and a typical 
such finite element function. 




Figure 2: Element diagrams. First row: discontinuous elements 
of degrees 0, 1, and 2 in two dimensions. Second row: Lagrange 
elements of degrees 1, 2, and 3 in two dimensions. Third and 
fourth rows: the corresponding elements in three dimensions. 



case, the face elements take as shape functions polynomial vector fields of the form 
p{x) = a + bx where a e M^, b gM. and x = {xi,X2,X3), a 4-dimensional subspace 
of the 12-dimensional space Pi(T, K^) of polynomial vector fields of degree at most 
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1. The degrees of freedom are taken to be the 0th order moments of the normal 
components on the faces of codimension 1, p Jjp{x) ■ Ufdx where / is a face 
and Hf the unit normal to the face. The element diagram is shown in the middle 
column of Figure ||. In the lowest order case the edge elements shape functions are 
polynomial vector fields of the form p{x) = a + b x x where a, 6 € M"^, which form 
a 6-dimensional subspace of Pi(r, M'^). The degrees of freedom are the 0th order 
moments over the edges of the component tangent to the edge, p i— > J^p{x) ■ dx, 
as indicated on the right of Figure 0. 




Figure 3: Element diagrams for some finite element approxima- 
tions to vector fields in two and three dimensions. Multiple dots 
are used as markers to indicate the evaluation of all components 
of a vector field. Arrows are used for normal moments on codi- 
mension 1 subsimplices and for tangential components on edges. 
Left: continuous piecewise linear fields. Middle: face elements of 
lowest order. Right: edge elements of lowest order. 

Each of these spaces can be generalized to arbitrarily high order. For the next 
higher order face space, the shape functions take the form p{x) — a{x) + b(x)x 
where a G Pi(T, R'^) and b S Pi(T) a linear scalar- valued polynomial. This gives a 
subspace oiF2{T,M.'^) of dimension 15, and the degrees of freedom are the moments 
of degree at most 1 of the normal components on the faces and the moments of 
degree of all components on the tetrahedron. This element is indicated on the left 
of Figure ^. For the second lowest order edge space, the shape functions take the 
form p{x) = a{x) + b{x) x x with a, 6 G Pi(T, M'^), giving a 20-dimensional space. 
The degrees of freedom are the tangential moments of degree at most 1 on the edges 
(two per edge) and the tangential moments of degree on the faces (two per face). 
This element is indicated on the right of Figure ^. 

The choice of the shape functions and the degrees of freedom determine the 
smoothness of the functions belonging to the assembled finite element space. For ex- 



142 



Douglas N. Arnold 




A 




Figure 4: The face (left) and edge (right) elements of the second 
lowest order in 2- and 3-dimensions. 



ample, the Lagrange finite element spaces of any degree belong to the Sobolev space 
H^{U) of L'^{fl) fimctions whose distributial first partial derivatives also belong to 
L^(f2) (and even to Unftyiyi)). In fact, the distributional first partial derivative of 
a continuous piecewise smooth function coincides with its derivative taken piecewise 
and so belongs to L^. Thus the degrees of freedom we imposed in constructing the 
Lagrange finite elements are sufficient to insure that the assembled finite element 
space Wh C if^(r2). In fact more is true: for the Lagrange finite element space with 
shape function spaces Wt = Tp{T), we have 

Wh — {u^ H^{Q) I u\t G Wt for all simplices T of the triangulation } . 

This says that, in a sense, the degrees of freedom impose exactly the continuity 
required to belong to H^, no less and no more. 

In contrast, the discontinuous piecewise polynomial spaces are subsets of L^(Q) 
but not of H^{^1), since their distributional first derivatives involve distributions 
supported on the interelement boundaries, and so do not belong to L^(f2). 

For the vector-valued finite elements there are more possibilities. The face and 
edge spaces contain discontinuous functions, and so are not contained in H^{^1, M.^). 
However, for vector fields belonging to one of the face spaces the normal component 
of the vector field does not jump across interelement boundaries, and this implies, 
via integration by parts, that the distributional divergence of the function coincides 
with the divergence taken piecewise. Thus the face spaces belong to H{div, fl), the 
space of vector fields on fl whose divergence belongs to L^. Indeed, for these 
spaces the degrees of freedom impose exactly the continuity of H{div), no less or 
more. For the edge spaces it can be shown that the tangential components of a 
vector field do not jump across element boundaries, and this implies that the edge 
functions belong to i?(curl, fi), the space of vector fields whose curl belongs to 
L^. Again the degrees of freedom impose exactly the continuity needed for inclusion 
in if (curl). 

3. Discrete differential complexes 

The de Rham complex 



d 



d 



d 



A"(0) - 
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is defined for an arbitrary smooth n-manifold fl. Here A "(^) denotes the space of 
differential A:-forms on il, i.e., for w e /\ (il) and x ^ il, uj{x) is an alternating k- 
linear map on the tangent space T^. The operators d : /\''~^^{^) denote 

exterior differentiation. This is is a complex in that the composition of two exterior 
differentiations always vanishes. Moreover, and if the manifold is topologically 
trivial, then it is exact. 

If is a domain in K'^, then we may identify its tangent space at any point 
with R^. Using the Euclidean inner product, the space of linear maps on R'^ may 
be identified by M.^ as usual, so A^(^) be identified with the space C°°(ri,IR^) 
of smooth vector fields on fl. Moreover, the space of alternating bilinear maps on 
may be identified with M'^ by associating to a vector u the alternating bilinear 
map {v,w) 1-^ det{u\v\'w). Thus we have an identification of A^(^) with M.^ as 
well. Finally the only alternating trilinear maps on are given by multiples of 
the determinant map i— > cdet(w|u|w), and so we may identify A'^(^) with 
C°°(f2). In terms of such proxy fields, the de Rham complex becomes 

C°°{n) C°°(f7,R3) C°°(f7,R3) '^'^ ! C°°(f^,R) ^ 0. 

(3.1) 

Alternatively we may consider i^-based forms and the sequence becomes 

H\n) Hicnvin) -^^^ H{diy,n) L^{n,R)~^0. 

The finite element spaces constructed above allow us to form discrete analogues 
of the de Rham complex. Given some triangulation of C R^, let Wh denote 
the space of continuous piecewise linear finite elements, Qh the lowest order edge 
element space, Sh the lowest order face element space, and Vh the space of piecewise 
constants. Then gradVF;i C Qh (since Qh contains all piecewise constant vector 
fields belonging to 7f(curl) and the gradient of a continuous piecewise linear is 
certainly such a function), cml Qh C Sh (since Sh contains all piecewise constant 
vector fields belonging to ff(curl)), and div^^i C Vh- Thus we have the discrete 
differential complex 

M Wh Qh Sh Vh ^ 0. (3.2) 

This differential complex captures the topology of the domain to the same extent 
as the de Rham complex. In particular, if the domain is topologically trivial, then 
the sequence is exact. 

It is convenient to abbreviate the above statement using the element diagrams 
introduced earlier. Thus we will say that the following complex is exact: 




By this we mean that if we assemble finite element spaces Wh , Qh , Sh , and Vh using 
the indicated finite elements and a triangulation of a topologically trivial domain, 
then the corresponding discrete differential complex ( |3.2[ ) is exact. 
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There is another important relationship between the de Rham complex (3.1) 
and the discrete complex (|]^). The defining degrees of freedom determine pro- 



W ; 



is 



jections : C°°(f}) ^ Wh, H)^ : C^i^.W) Qh, and so on. In fact 
just the usual interpolant, U]^ is the L^-projection into the piecewise constants, 
and the projections 11^ and onto the edge and face elements are determined by 
the maintenance of the appropriate moments. It can be checked, based on Stokes 
theorem, that the following diagram commutes. 



grad 



Qh 



curl 



curl 



Sh 



div 



div 



^ 
(3.3) 

The finite element spaces appearing in this diagram, with one degree of freedom 
for each vertex for Wh, for each edge for Qh, for each face for Sh, and for each 
simplex for Vh, are highly geometrical. In fact, recalling the identifications between 
fields and differential forms, we may view these spaces as spaces of piecewise smooth 
differential forms. They were in fact first constructed in this context, without any 
thought of finite elements or numerical methods, by Whitney |^3|. The spaces were 
reinvented, one-by-one, as finite element spaces in response to the needs of various 
numerical problems, and the properties which are summarized in the commutative 
diagram above were slowly rediscovered as needed to analyze the resulting numerical 
methods. The connection between low order edge and face finite elements and 
Whitney forms was first realized by Bossavit . 

Analogous statements hold for higher order Lagrange, edge, face, and discon- 
tinuous finite elements. For example, the following diagram commutes and has 
exact rows: 



rad 








lad 






We shall see many other discrete differential complexes below. 



4. Stability of Galerkin methods 

Consider first the solution of the Dirichlet problem for Poisson's equation on 
a domain in R": 

—Au — f in n, u = on d^l. 
The solution can be characterized as the minimizer of the energy functional 

1 



£iu) 



|gradM(x)| dx — f{x)u{x)dx 
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over the Sobolev space H^{il) (consisting of H^{il) fmictions vanishing on dil), or 
as the solution of the weak problem: find u e H^{il.) such that 

/ grad.(x) . grad.(x) = [ for all . e H\n). 

Jn Jn 

We may define an approximate solution Uh by minimizing the Dirichlet integral 
over a finite dimensional subspace Wh of H^{^1); this is the classical Ritz method. 
Equivalently, we may use the Galerkin method, in which G Wh is determined by 
the equations 



gradM;i(a;) • gradw(a:) dx = f{x)v{x) dx for all v £ Wh- 
n Jn 

After choice of a basis in Wh this leads to a system of linear algebraic equations, 
and Uh is computable. 

Let Th denote the discrete solution operator f Uh- Then it is easy to check 
that Th is bounded as a linear operator from H^^{^l) := H^{il)* to H^{^1) by a 
constant that depends only on the domain (and, in particular, doesn't increase 
if the space Wh is enriched). This says that the Galerkin method is stable. A 
consequence is the quasioptimality estimate 

\\u-Uh\\m <c inf \\u-v\\hi, (4.4) 
veWh 

for some constant c depending only on the domain fl. Note that there is no restric- 
tion on the subspace Wh to obtain this estimate. Galerkin's method for a coercive 
elliptic problem is always stable and convergence depends only on the approxima- 
tion properties of the subspace. A natural choice for Wh is the Lagrange finite 
element space of some degree p with respect to some regular simplicial mesh of 
maximal element size /i, in which case Galerkin's method is a standard finite ele- 



ment method. In this case the right hand side of (|4.4| ) is 0{hP) provided that u is 
sufficiently smooth. 

Next consider the related eigenvalue problem, which arises in the determination 
of the fundamental frequencies of a drum. That is, we seek standing wave solutions 
w{x,t) to the wave equation on some bounded domain f2 C which vanish on 
dfl. Assuming that the tension and density of the drum membrane are unity, these 
solutions have the form w{x, t) = a coa{\/\t)u{x) + (3 s\n{-/\t)u{x) where a and /3 
are constants and u and A satisfy the eigenvalue problem 

— Au — \u \n ^l, u = on dVL. 

The eigenvalues A form a sequence of positive numbers tending to infinity. The 
numbers a/A/ (27r) are the fundamental frequencies of the drum and the functions 
u give the corresponding fundamental modes. 

The eigenvalues and eigenfunctions are characterized variationally as the crit- 
ical values and critical points of the Rayleigh quotient 

_ S^\gra.du{x)\^ dx 
'^^"^ S^Hx)\^dx ' 
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defined for nonzero u belonging to the Sobolev space H^{^). The classical Rayleigh- 
Ritz method for the approximation of eigenvalue problems determines approximate 
eigenvalues Xh and eigenfunctions Uh as the critical values and points of the restric- 
tion of TZ to the nonzero elements of some finite dimensional subspace Wt of H^(yi). 
Equivalently, we can write the eigenvalue problem in weak form: find A G M and 
nonzero u G if^(ri) such that 

rad.(.) . grad.(x) d. = X ( d. for all . G H\n). (4.5) 

The Galerkin approximation of the eigenvalue problem, which is equivalent to the 
Rayleigh-Ritz method, seeks A/i G M and nonzero Uh G Wh such that 



Taduh{x) ■ gradw(a;) dx — Xh / Uh{x)v{x) dx for all v G Wh- (4-6) 
vt Jn 

We now discuss the convergence of this method. Let A denote the jth eigen- 



value of the problem (4.5). In the interest of simplicity we assume that A is a simple 
eigenvalue, so the corresponding eigenfunction u is uniquely determined up to sign 
by the normalization — 1- Similarly let A^ and denote the jth eigenvalue 



of (4.6). It can then be proved (see, e.g., H for much more general results) that 



there exists a constant c such that 

11" - w/illfl-i < c inf \\u~v\\hi, |A - A;i| < c||u - (4.7) 
veWh 

In short, the eigenfunction approximation is quasioptimal and the eigenvalue error 
is bounded by the square. Again there is no restriction on the space Wh- 

Figure |^ reports on the computation of the eigenvalues of the Laplacian on an 
elliptical domain of aspect ratio 3 using Lagrange finite elements of degree 1. 

Now consider an analogous problem, the computation of the resonant frequen- 
cies of an electromagnetic cavity occupying a region f2 C M'^. In this case we wish 
to find standing wave solutions of Maxwell's equations. If we take the electric per- 
mittivity and the magnetic permeability to be unity and assume a lossless cavity 
with perfectly conducting boundary, we are led to the following eigenvalue problem 
for the electric field: find nonzero E : il ^ M.^, A G K. such that 

curl curl £: = A£', div£; = Oinf}, E x n = on dn- (4.8) 

This is again an elliptic eigenvalue problem and the eigenvalues form a sequence of 
positive numbers tending to infinity. The divergence constraint is nearly redundant 
in this eigenvalue problem. Indeed if curl curl i? = XE for A > 0, then div E = 
X^^ div curl curl i? = since the divergence of a curl vanishes. Thus the eigenvalue 
problem 

cmlcmlE = XE iiifl, E x n ^ on dn, (4.9) 



has the same eigenvalues and eigenfunctions as ( [4.8| ) except that it also admits 
A = as an eigenvalue, and the corresponding eigenspace is infinite-dimensional (it 
contains the gradients of all smooth functions vanishing on the boundary of fi). The 
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Figure 5: The point plot shows the first 40 eigenvalues computed 
with piecewise linear finite elements with respect to the triangu- 
lation shown (•) versus the exact eigenvalues (+). The surface 
plot shows the computed eigenfunction associated to the fourth 
eigenvalue. The mesh has 737 vertices, of which 641 are interior, 
and 1,376 triangles. 



eigenvalues and eigenfunctions are now critical points and values of the Rayleigh 
quotient 



n{E) 



^^\cnv\E{x)\'^dx 



over the space of nonzero fields E in i7(curl, 51), which is defined to be the space 
of functions for which both the above integrals exist and are finite and which have 
vanishing tangential component on the boundary (i.e., i5 x n = on d^). 

In Figure ^ we show the result of approximating a two-dimensional version 
of this eigenvalue problem using the Rayleigh-Ritz method or, equivalently, the 
Galerkin method with continuous piecewise linear vector fields on whose tangen- 
tial components vanish on the boundary (the first element depicted in Figure 
For r2 we take a square of side length tt, in which case the nonzero eigenvalues are 
known to be all numbers of the form A = rn? + with < m^n G 7^ not both 
zero, and the corresponding eigenfunctions are E — (sin my, sin nx). For the mesh 
pictured, the finite element space has dimension 290. We find that 73 of the 290 
computed eigenvalues are between and 10 and that they have no tendency to clus- 
ter near the integers 1, 1, 2, 4, 4, 5, 5, 8, 9, 9 which are the exact eigenvalues between 
and 10. Thus this numerical method is useless: th e c omputed eigenvalues bear 
no relation to the true eigenvalues! The analogue of ( [l.7| ) is surely not true. 

If instead we choose the lowest order edge elements as the finite element space 
(Figure ^, top right), we get very different results. Using the same mesh, the edge 
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Figure 6: The plot shows the first 73 eigenvalues computed with 
piecewise linear finite elements for the resonant cavity problem on 
the square using the mesh shown. They bear no relation to the 
exact eigenvalues, 1, 1, 2, 4, 4, ... , indicated by the horizontal 
lines. 



finite element space has dimension 472. It turns out that 145 of the computed eigen- 
values are zero (to within round-off), and the subsequent eigenvalues are 0.9998, 
0.9999, 2.0023, 3.9968, 4.0013, . . . , i.e., excellent approximations of the exact eigen- 
values. See Figure 0. 




Figure 7: The first plot shows the first 100 positive eigenvalues for 
the resonant cavity problem on the square computed with lowest 
order edge elements using the mesh of Figure ||. The error m 
the first 54 eigenvalues is below 2%. The inset focuses on the 
first 10 eigenvalues, for which the error is less than 0.25%. The 
second plot shows the vector field associated to the third positive 
eigenvalue. 
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The striking difference between the behavior of the continuous piecewise linear 
finite elements and the edge elements for the resonant cavity problem is a question 
of stability. We shall return to this below, after examining stability in a simpler 
context. 



5. Stability of mixed formulations 

Consider now the Dirichlet problem 

— divC gradw — f in fl, u = on 957, 

where $7 is a domain in R'^ and the coefficient C is a symmetric positive definite 
matrix at each point. We may again characterize m as a minimizer of the energy 
functional 

u t-^ — C grad u ■ grad udx — / fudx 



and use the Ritz method. This procedure is always stable. 

However, for some purposes it is preferable to work with the equivalent first 
order system 

CT = Cgradu, — diver = /. (5.10) 

The pair (ct, u) is then characterized variationally as the unique critical point of the 
functional 

C{a,u) = / {—C~^(T-(T + udiva)dx'' / fudx (5-11) 
Jn 2 Jo 

over H{div,n) x L^(f7). Note that {a,u) is a saddle-point of C, not an extremum. 
Numerical discretizations based on such saddle-point variational principles are called 
mixed methods. 



It is worth interpreting the system (5.10) in the language of differential forms, 



because this brings some insight. The function u is a 0-form, and the operation 
u I— *■ gradu is just exterior differentiation. The vector field ct is a proxy for a 2-form 
and the operation a i— *■ div a is again exterior differentiation. The loading function 
/ is the proxy for a 3-form. Since gradu is the proxy for a 1-form, it must be 
that the operation on differential forms that corresponds to multiplication by C 
takes 1-forms to 2-forms. In fact, if we untangle the identifications, we find that 
multiplication by C is a Hodge star operation. A Hodge star operator defines an 
isomorphism of onto To determine a particular such operator, we 

must define an inner product on the tangent space at each point of fl. The pos- 
itive definite matrix C does exactly that. Many of the partial differential equations 
of mathematical physics admit similar interpretations in terms of differential forms. 
For a discussion of this in the context of discretization, see . 

A natural approach to discretization of the mixed variational principle is to 
choose subspaces Sh C i?(div, 17), Vh C i^(ri) and seek a critical point {ah,Uh) & 
Sfi^Vfi- This is of course equivalent to a Galerkin method and leads to a system of 
linear algebraic equations. However in this case, stability is not automatic. It can 
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happen that the discrete system is singular, or more commonly, that the norm of 
the discrete solution operator grows unboundedly as the mesh is refined. 

In a fundamental paper, Brezzi |^ established two conditions that together 
are sufficient (and essentially necessary) for stability. Brezzi's theorem applied to 
a wide class of saddle-point problems, but for simplicity we will state the stability 



conditions for the saddle-point problem associated to the functional ( 5.11 ) 
(SI) There exists 71 > such that 



/ 



C V -Tdx > 7l|k|lH(div)' 



for all r S Sh such that / div tv dx ~ for all v E Vh- 
(S2) There exists 72 > such that for all v E Vh there exists nonzero t E Sh 
satisfying 



Jn 



vdiyrdx > 72|lw|lL2|lT||H(dr 



Theorem (Brezzi) // the stability conditions (SI) and (S2) are satisfied, then C 
admits a unique critical point (ah,Uh) over Sh x Vh, the solution operator f i—f 
{ah,Uh) is bounded i^(il) — > H{drv,Q) x L^(ri), and the quasioptimal estimate 

Ik - cr/i||H(div) + < c inf (||ct - T||H(div) + ||m - w||l2) 

holds with c depending on 71 and 72. 

The stability conditions of Brezzi strongly limit the choice of the mixed finite 
element spaces Sh and Vh- Condition (SI) is satisfied if the indicated functions 
T E Sh, those whose divergence is orthogonal to Vh, are in fact divergence- free. (In 
practice, this is nearly the only way it is satisfied.) This certainly holds if div 5/1 C 
Vh, and so such as inclusion is a common design principle of mixed finite element 
spaces. On the other hand, condition (S2) is most easily satisfied if div5/i D Vh, 
because in this case, given v E Vh, we can choose r E Sh with divr = v, so 
V div rdx — ||i'|j^2, and the second condition will be satisfied as long as we can 
insure that ||T||//(div) < 7;7^II''^IU2- short, we need to know that div maps Sh 
onto Vh and that div jg^ admits a bounded one-sided inverse. 

The face elements of Raviart-Thomas and Nedelec were designed to satisfy 
both these conditions. Specifically, let Sh again denote the space of face elements of 
lowest degree (whose element diagram is shown in the middle of the second row of 
Figure and Vh the space of piecewise constants.]^ We know that Sh C if (div, ft) 
so these elements are admissable for the mixed variational principle. Moreover, we 
have div Sh C Vh, so (SI) holds. 

To verify (S2), we refer to the commutative diagram (^^). Given v E Vh, we 
can solve the Poisson equation A0 = v and take a — grad cj) to obtain a function 



^It may seem odd to seek Uf^ in Vji, a space of discrete 3-forms, rather than in a space of 
0-forms, since u is a 0-form. The resolution is through a Hodge star operator, this time formed 
with respect to the EucUdean inner product on R^. In the mixed method u/^ is a discrete 3-form, 
approximating the image of u under this star operator. 
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with divfT — V and HctHhi ^ C\\v\\l2. Now let t — n^cr G Sh- Then 



div T = div n^(T = div a ~ Ilf^ v ^ v, 
where we have used the commutativity and the fact that v e Vh- Moreover 



Lff(div; 



< c||cr||/fi < c'||u||^2, where we used the boundedness of onH^{Q, 



This shows that div Vh — Sh and establishes a bound on the one-sided inverse, and 
so verifies (S2). Of course, the same argument shows the stability of a mixed method 
based on higher order face elements as well. 

Thus we see that the stability of the mixed finite element method depends 
on the properties of the spaces Vh and Sh encoded in the rightmost square of the 
commutative diagram (O). 



Now let us return to the resonant cavity eigenvalue problem (4.9) for which 
we explored the Galerkin method: find A/i e M, 7^ G Qh such that 



curl_E;i • curlFdx — Xh / Eh ■ F dx for all F e Qh 



(5.12) 



We saw that it Qh C 7? (curl, il) is taken to be a space of edge elements this method 
gives good results in that the positive eigenvalues of the discrete problem are good 
approximations for the positive eigenvalues of the continuous problem. However, 
the simple choice of Lagrange finite elements did not give good results. We now 
explain the good performance of the edge elements based on the middle square 
of the commutative diagram (3.3). Following Boffi et. al Q we set Ph = curl Qh 
and introduce the following mixed discrete eigenvalue problem: find A/i G M, 7^ 
{Eh,Ph) & Qh^ Pfi such that 



Eh-Fdx 



curl F ■ Ph dx — for all F £ Qh 



(5.13) 



cmlEh ■ qdx = —Xh / Ph ■ qdx for all q G Ph- (5-14) 
n Jn 

It is then easy to verify that if Xh, Eh is a solution to (^.12[ ) with Xh > 0, then 



Xhj{Eh, Xf^^ cml Eh) is a solution to (5.13), and if Xh, {Eh,Ph) is a solution to 



( 5.13| ) then Xh > and Xh, Eh is a solution to (5.12). In short, the two problems 
are equivalent except that the former admits a zero eigenspace which the mixed 
formulation suppresses. As explained in Q, the accuracy of the mixed eigenvalue 
problem ( 5.13| ) hinges on the stability of the corresponding mixed source problem. 
This is a saddle-point problem of the sort studied by Brezzi, and so stability depends 
on conditions analogous to (SI) and (S2). The proof of these conditions in case Qh 
is the space of edge elements follows, as in the preceding stability verification, from 



surjectivity and commutativity properties encoded in the diagram (3.3). 

The diagram can also be used to explain the zero eigenspace computed with 
edge elements. Recall that in the case of the mesh shown in Figure this space had 
dimension 145. In fact, this eigenspace is simply the null space of t he c url operator 
restricted to Qh- Referring again to the commutative diagram (3^), this is the 
gradient of the space Wh of linear Lagrange elements vanishing on the boundary. 
Its dimension is therefore exactly the number of interior nodes of the mesh. 
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6. The elasticity complex 

Let § denote the space of 3 x 3 symmetric matrices. Given a volumetric loading 
density f : fl ^ M.^, the system of linearized elasticity determines the displacement 
field w : O ^ and the stress field cr : ^ S induced in the elastic domain fl by 
the equations 

a ~ C eu, — div a = f, 

together with boundary conditions such as m = on 90. Here eu is the symmetric 
part of the matrix gradu, and the elasticity tensor C : § ^ § is a symmetric positive 
definite linear operator describing the particular elastic material, possibly varying 
from point to point. 

The solution (ct, u) may be characterized variationally as a saddle-point of the 
Hellinger-Reissner functional 

C{(T,u) — j { — C^^ (7 : a + u ■ div a)dx ~ j f-udx (6.15) 

over iJ(div, 51,§) x i^(0,R^) (i.e., a is sought in the space of square-integrable 
symmetric-matrix-valued functions whose divergence by rows is square-integrable, 
and u is sought among all square-integrable vector fields). 

For a mixed finite element method, we need to specify finite element subspaces 
Sh C i?(div, 0,§) and Vh C L^(0,]R^) and restrict the domain of the variational 
problem. Of course the spaces must be carefully designed if the mixed method 
is to be stable: the analogues of the stability conditions (SI) and (S2) must be 
satisfied. The functional ( |6.15D is quite similar in appearance to (5.11) and so 



it might be expected that the mixed finite elements developed for the latter (the 
face elements for a and discontinuous elements for u) could be adapted to the 
case of elasticity. In fact, the requirement of symmetry of the stress tensor and, 
correspondingly, the replacement of the gradient by the symmetric gradient, changes 
the structure significantly. Four decades of searching for mixed finite elements for 
elasticity beginning in the 1960s did not yield any stable elements with polynomial 
shape functions. 

Using discrete differential complexes, R. Winther and the author recently de- 
veloped the first such elements for elasticity problems in two dimensions ||l[. (The 
three-dimensional case remains open.) For elasticity, the displacement and stress 
fields cannot be naturally interpreted as differential forms and the relevant differ- 
ential complex is not the de Rham complex. In three dimensions it is instead the 
elasticity complex: 

T-^C°°(0,M3) — C°°(0,§) — ^ C°°(0,§) """^ ) C°°(0,R3)^0. 

Here the operator J is a second order differential operator which acts on a symmetric 
matrix field by first replacing each row with its curl and then replacing each column 
with its curl to obtain another symmetric matrix field. The resolved space T is the 
six-dimensional space of infinitesimal rigid motions, i.e., the same space of linear 
polynomials a + b x x which arose as the shape functions for the lowest order edge 
elements. If the domain is topologically trivial, this complex is exact. Although 
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it involves a second order differential operator, and so looks quite different from the 
de Rham complex, Eastwood recently pointed out that it can be derived from 
the de Rham complex via a general construction known as the Bernstein-Gelfand- 
Gelfand resolution. 

In two dimensions the elasticity complex takes the form 



Pi ^c°°{n) 



C°°(f7,§) 



0, 



where now the second order differential operator is 



J = 



dxl 



92 



\ dxidx2 



92 \ 



dxidx2 



dxl 




Figure 8: Element diagram for the new mixed finite elements for 
elasticity, lowest order case. 



In the lowest order case, the finite elements we introduced in for which 
the element diagrams can be seen in Figure |8[ use discontinuous piecewise linear 
vector fields for the displacement field and a piecewise polynomial space which we 
shall now describe for the stress field. The shape functions on an arbitrary triangle 
T are given by 

5T = {reP3(T,§)| divrePi(r,M2)}^ 

which is a 24-dimensional space consisting of all quadratic symmetric matrix fields 
on T together with the divergence-free cubic fields. The degrees of freedom are 

• the values of three components of t{x) at each vertex a; of T (9 degrees of 
freedom) 

• the values of the moments of degree and 1 of the two components of rn on 
each edge e of T (12 degrees of freedom) 

• the value of the three components of the moment of degree of r on T (3 
degrees of freedom) 

Note that these degrees of freedom are enough to ensure continuity of rn across 
element faces, and so will furnish a finite element subspace of iJ(div, i7, §). The 
continuity is not however, the minimal needed for inclusion in _ff(div). The de- 
grees of freedom also enforce continuity at the vertices, which is not required for 
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membership in _ff(div). For various reasons, it would be useful to have a mixed 
finite element for elasticity that does not use vertex degrees of freedom. But, as we 
remark below, this is not possible if we restrict to polynomial shape functions. 

In order to have a well-defined finite element, we must verify that the 24 degrees 
of freedom form a basis for the dual space of St- We include this verification since 
it illustrates an aspect of the role of the elasticity complex. Since dim St = 24, 
we need only show that if all the degrees of freedom vanish for some t £ St, then 
r = 0. Now rn varies cubically along each edge, vanishes at the endpoints, and has 
vanishing moments of degree and 1. Therefore rn = 0. Letting v = divr, a linear 
vector field on T, we get by integration by parts that 

dx = — I T : ev dx + / rn ■ v ds ~ 
JT Jar 

since the integral of r vanishes as well as rn. Thus r is divergence-free. In view of 
the exactness of the elasticity complex, t — Jq for some smooth function q. Since 
all the second partial derivatives of q belong to P3(r), q S F^lT). Adjusting by 
an element of Pi(r) (the null space of J), we may take q to vanish at the vertices. 
Now d^q/ds^ — rn ■ n = Q on each edge, whence q is identically zero on dT. This 
implies that the gradient of q vanishes at the vertices. Since d'^q/dsdn = —rn-t = 
on each edge (with t a unit vector tangent to the edge), we conclude that dq/dn 
vanishes identically on dT as well. Since q has degree at most 5, it must vanish 
identically. 

Let : C°°(ri,§) Sh denote the projection associated with the supplied 
degrees of freedom, and Ii\ : C°°(il,M^) Vh the L^-projection. For any triangle 
T, r e C°°(r2,§), and v G Pi(f7,M2)^ ^e have 

/ div(r - nfr) ■ v dx = - I (r - II^t) : evdx+ (r - Ilf r)n • vds. 

JT JT JOT 

The degrees of freedom entering the definition of 11^ ensure that the right hand side 
vanishes, and from this we obtain the commutativity divll^r = 11^ divr which 
is essential for stability. (Actually a technical difficulty arises here, since 11^ as 
given is not bounded on _ff^(r2,§). See for the resolution.) Note that, by their 
definitions, div5/i C Vh and, using the commutativity, we have div^^, = Vh, i.e., 
Sh — "^'^ > V/i — > is exact. To complete this to a discrete analogue of the elasticity 
complex, we define Yh to be the inverse image of Sh under J. Then Yh is exactly 
the space of piecewise quintic polynomials which are at the vertices of the 
meshes. This is in fact a well-known finite element space, called the Hermite quintic 
or Argyris space, developed for solving 4th order partial differential equations (for 
which the inclusion in H^{fl) and therefore continuity is required). The shape 
functions are P5(T) and the 21 degrees of freedom are the values of the function 
and all its first and second partial derivatives at the vertices and the integrals of 
the normal derivatives along edges. We then have a discrete elasticity complex 

Pi ^Yh — ^ Sh Vh^ 0, 
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or, diagrammatically, 





Moreover this sequence is exact and is coupled to the two-dimensional elasticity 
sequence via a commuting diagram: 



Pi ^c°^{n) 



C°°(fl,§) 



Vh - 











The right half of this diagram encodes the information necessary to establish the 
stability of our mixed finite element method. 

The Hermite quintic finite elements arose naturally from our mixed finite ele- 
ments to complete the commutative diagram. Had they not been long known, we 
could have used this procedure to devise a finite element space contained in H^{n). 
In fact, on close scrutiny we can see that any stable mixed finite elements for elas- 
ticity with polynomial shape functions will give rise to a finite element space with 
polynomial shape functions contained in H^{il). However, it is known that such 
spaces are difficult to construct and complicated. In fact, it can be proved that an 
finite element space must utilize shape functions of degree at least 5 and the first 
and second partial derivatives at the vertices must be among the degrees of freedom 
jl^ . This helps explain why mixed finite elements for elasticity have proven so hard 
to devise. In particular, we can rigorously establish the stress elements must involve 
polynomials of degree 3, and that vertex degrees of freedom are unavoidable. 

In addition to the element just described, elements of all greater orders are 
also introduced in The elements of next higher order can be seen as the final 
two elements in this discrete elasticity complex. 





It is also possible to simplify the lowest order element slightly. To do this we reduce 
the displacement space from piecewise linear vector fields to piecewise rigid motions, 
and we replace the stress space with the inverse image under the divergence of the 
reduced displacement space. This leads to a stable element shown in this exact 
sequence: 





Because of the unavoidable complexity of finite elements, practitioners 
solving 4th order equations often resort to nonconforming finite element approxi- 
mations of H^. This means that the finite element space does not belong to 
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in that the function or the normal derivative may jump across element boundaries, 
but the spaces are designed so that jumps are small enough in some sense (e.g., on 
average). The error analysis is more complicated for nonconforming elements, since 
in addition to stability and approximation properties of the finite element space, one 
must analyze the consistency error arising from the jumps in the finite elements. 
In 1^ Winther and the author investigated the the possibility of nonconforming 
mixed finite elements for elasticity, which, however are stable and convergent, and 
developed two such elements. These are related to nonconforming H'^ elements via 
nonconforming discrete elasticity complexes, two of which arc pictured here: 




In both cases the shape function space for the stress is contained between Pi (T, §) 
and P2 (T, S) . The nonconforming finite element depicted in these diagrams 
was developed for certain 4th order problems in Note the nonconforming 

mixed elasticity elements are significantly simpler than the conforming ones (and, 
in particular, don't require vertex degrees of freedom). 
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